{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "@import url('http://fonts.googleapis.com/css?family=Source+Code+Pro');\n",
       "@import url('http://fonts.googleapis.com/css?family=Vollkorn');\n",
       "@import url('http://fonts.googleapis.com/css?family=Arimo');\n",
       "@import url('http://fonts.googleapis.com/css?family=Fira_sans');\n",
       "\n",
       "    div.cell{\n",
       "        width: 900px;\n",
       "        margin-left: 0% !important;\n",
       "        margin-right: auto;\n",
       "    }\n",
       "    div.text_cell code {\n",
       "        background: transparent;\n",
       "        color: #000000;\n",
       "        font-weight: 600;\n",
       "        font-size: 11pt;\n",
       "        font-style: bold;\n",
       "        font-family:  'Source Code Pro', Consolas, monocco, monospace;\n",
       "   }\n",
       "    h1 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "\t}\n",
       "\t\n",
       "    div.input_area {\n",
       "        background: #F6F6F9;\n",
       "        border: 1px solid #586e75;\n",
       "    }\n",
       "\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 30pt;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h2 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "        text-align: left;\n",
       "    }\n",
       "    .text_cell_render h2 {\n",
       "        font-weight: 200;\n",
       "        font-size: 16pt;\n",
       "        font-style: italic;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1.5em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h3 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h3 {\n",
       "        font-weight: 200;\n",
       "        font-size: 14pt;\n",
       "        line-height: 100%;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 2em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    }\n",
       "    h4 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h4 {\n",
       "        font-weight: 100;\n",
       "        font-size: 14pt;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    h5 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 200;\n",
       "        font-style: normal;\n",
       "        color: #1d3b84;\n",
       "        font-size: 16pt;\n",
       "        margin-bottom: 0em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    div.text_cell_render{\n",
       "        font-family: 'Fira sans', verdana,arial,sans-serif;\n",
       "        line-height: 125%;\n",
       "        font-size: 115%;\n",
       "        text-align:justify;\n",
       "        text-justify:inter-word;\n",
       "    }\n",
       "    div.output_subarea.output_text.output_pyout {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_subarea.output_stream.output_stdout.output_text {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_wrapper{\n",
       "        margin-top:0.2em;\n",
       "        margin-bottom:0.2em;\n",
       "}\n",
       "\n",
       "    code{\n",
       "      font-size: 70%;\n",
       "    }\n",
       "    .rendered_html code{\n",
       "    background-color: transparent;\n",
       "    }\n",
       "    ul{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li li{\n",
       "        padding-left: 0.2em; \n",
       "        margin-bottom: 0.2em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    ol{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ol li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    a:link{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:visited{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:hover{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:focus{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:active{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    .rendered_html :link {\n",
       "       text-decoration: underline; \n",
       "    }\n",
       "    .rendered_html :hover {\n",
       "       text-decoration: none; \n",
       "    }\n",
       "    .rendered_html :visited {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :focus {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :active {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "    } \n",
       "    hr {\n",
       "      color: #f3f3f3;\n",
       "      background-color: #f3f3f3;\n",
       "      height: 1px;\n",
       "    }\n",
       "    blockquote{\n",
       "      display:block;\n",
       "      background: #fcfcfc;\n",
       "      border-left: 5px solid #c76c0c;\n",
       "      font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "      width:680px;\n",
       "      padding: 10px 10px 10px 10px;\n",
       "      text-align:justify;\n",
       "      text-justify:inter-word;\n",
       "      }\n",
       "      blockquote p {\n",
       "        margin-bottom: 0;\n",
       "        line-height: 125%;\n",
       "        font-size: 100%;\n",
       "      }\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    scale:100,\n",
       "                        availableFonts: [],\n",
       "                        preferredFont:null,\n",
       "                        webFont: \"TeX\",\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# format the book\n",
    "%matplotlib inline\n",
    "import sys\n",
    "from __future__ import division, print_function\n",
    "import sys\n",
    "sys.path.insert(0,'../code')\n",
    "import book_format\n",
    "book_format.load_style('../code')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Observer Bias\n",
    "\n",
    "## The Red Line problem\n",
    "\n",
    "In Massachusetts, the Red Line is a subway that connects Cambridge and\n",
    "Boston. When I was working in Cambridge I took the Red Line from Kendall\n",
    "Square to South Station and caught the commuter rail to Needham. During\n",
    "rush hour Red Line trains run every 7–8 minutes, on average.\n",
    "\n",
    "When I arrived at the station, I could estimate the time until the next\n",
    "train based on the number of passengers on the platform. If there were\n",
    "only a few people, I inferred that I just missed a train and expected to\n",
    "wait about 7 minutes. If there were more passengers, I expected the\n",
    "train to arrive sooner. But if there were a large number of passengers,\n",
    "I suspected that trains were not running on schedule, so I would go back\n",
    "to the street level and get a taxi.\n",
    "\n",
    "While I was waiting for trains, I thought about how Bayesian estimation\n",
    "could help predict my wait time and decide when I should give up and\n",
    "take a taxi. This chapter presents the analysis I came up with.\n",
    "\n",
    "This chapter is based on a project by Brendan Ritter and Kai Austin, who\n",
    "took a class with me at Olin College. The code in this chapter is\n",
    "available from <http://thinkbayes.com/redline.py>. The code I used to\n",
    "collect data is in <http://thinkbayes.com/redline_data.py>. For more\n",
    "information see Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The model\n",
    "\n",
    "![PMF of gaps between trains, based on collected data, smoothed by KDE.\n",
    "`z` is the actual distribution; `zb` is the biased\n",
    "distribution seen by passengers. ](figs/redline0.pdf)\n",
    "\n",
    "[fig.redline0]\n",
    "\n",
    "Before we get to the analysis, we have to make some modeling decisions.\n",
    "First, I will treat passenger arrivals as a Poisson process, which means\n",
    "I assume that passengers are equally likely to arrive at any time, and\n",
    "that they arrive at an unknown rate, $\\lambda$, measured in passengers\n",
    "per minute. Since I observe passengers during a short period of time,\n",
    "and at the same time every day, I assume that $\\lambda$ is constant.\n",
    "\n",
    "On the other hand, the arrival process for trains is not Poisson. Trains\n",
    "to Boston are supposed to leave from the end of the line (Alewife\n",
    "station) every 7–8 minutes during peak times, but by the time they get\n",
    "to Kendall Square, the time between trains varies between 3 and 12\n",
    "minutes.\n",
    "\n",
    "To gather data on the time between trains, I wrote a script that\n",
    "downloads real-time data from\n",
    "<http://www.mbta.com/rider_tools/developers/>, selects south-bound\n",
    "trains arriving at Kendall square, and records their arrival times in a\n",
    "database. I ran the script from 4pm to 6pm every weekday for 5 days, and\n",
    "recorded about 15 arrivals per day. Then I computed the time between\n",
    "consecutive arrivals; the distribution of these gaps is shown in\n",
    "Figure [fig.redline0], labeled `z`.\n",
    "\n",
    "If you stood on the platform from 4pm to 6pm and recorded the time\n",
    "between trains, this is the distribution you would see. But if you\n",
    "arrive at some random time (without regard to the train schedule) you\n",
    "would see a different distribution. The average time between trains, as\n",
    "seen by a random passenger, is substantially higher than the true\n",
    "average.\n",
    "\n",
    "Why? Because a passenger is more like to arrive during a large interval\n",
    "than a small one. Consider a simple example: suppose that the time\n",
    "between trains is either 5 minutes or 10 minutes with equal probability.\n",
    "In that case the average time between trains is 7.5 minutes.\n",
    "\n",
    "But a passenger is more likely to arrive during a 10 minute gap than a 5\n",
    "minute gap; in fact, twice as likely. If we surveyed arriving\n",
    "passengers, we would find that 2/3 of them arrived during a 10 minute\n",
    "gap, and only 1/3 during a 5 minute gap. So the average time between\n",
    "trains, as seen by an arriving passenger, is 8.33 minutes.\n",
    "\n",
    "This kind of **observer bias** appears in many contexts.\n",
    "Students think that classes are bigger than they are because more of\n",
    "them are in the big classes. Airline passengers think that planes are\n",
    "fuller than they are because more of them are on full flights.\n",
    "\n",
    "In each case, values from the actual distribution are oversampled in\n",
    "proportion to their value. In the Red Line example, a gap that is twice\n",
    "as big is twice as likely to be observed.\n",
    "\n",
    "So given the actual distribution of gaps, we can compute the\n",
    "distribution of gaps as seen by passengers. `BiasPmf` does\n",
    "this computation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def BiasPmf(pmf):\n",
    "    new_pmf = pmf.Copy()\n",
    "\n",
    "    for x, p in pmf.Items():\n",
    "        new_pmf.Mult(x, x)\n",
    "\n",
    "    new_pmf.Normalize()\n",
    "    return new_pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pmf` is the actual distribution; `new_pmf` is the biased\n",
    "distribution. Inside the loop, we multiply the probability of each\n",
    "value, `x`, by the likelihood it will be observed, which is\n",
    "proportional to `x`. Then we normalize the result.\n",
    "\n",
    "Figure [fig.redline0] shows the actual distribution of gaps, labeled\n",
    "`z`, and the distribution of gaps seen by passengers, labeled\n",
    "`zb` for “z biased”."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wait times\n",
    "\n",
    "![CDF of `z`, `zb`, and the wait time seen by\n",
    "passengers, `y`. ](figs/redline2.pdf)\n",
    "\n",
    "[fig.redline2]\n",
    "\n",
    "Wait time, which I call `y`, is the time between the arrival\n",
    "of a passenger and the next arrival of a train. Elapsed time, which I\n",
    "call `x`, is the time between the arrival of the previous\n",
    "train and the arrival of a passenger. I chose these definitions so that\n",
    "`zb = x + y`.\n",
    "\n",
    "Given the distribution of `zb`, we can compute the\n",
    "distribution of `y`. I’ll start with a simple case and then\n",
    "generalize. Suppose, as in the previous example, that `zb` is\n",
    "either 5 minutes with probability 1/3, or 10 minutes with probability\n",
    "2/3.\n",
    "\n",
    "If we arrive at a random time during a 5 minute gap, `y` is\n",
    "uniform from 0 to 5 minutes. If we arrive during a 10 minute gap,\n",
    "`y` is uniform from 0 to 10. So the overall distribution is a\n",
    "mixture of uniform distributions weighted according to the probability\n",
    "of each gap.\n",
    "\n",
    "The following function takes the distribution of `zb` and\n",
    "computes the distribution of `y`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def PmfOfWaitTime(pmf_zb):\n",
    "    metapmf = thinkbayes.Pmf()\n",
    "    for gap, prob in pmf_zb.Items():\n",
    "        uniform = MakeUniformPmf(0, gap)\n",
    "        metapmf.Set(uniform, prob)\n",
    "\n",
    "    pmf_y = thinkbayes.MakeMixture(metapmf)\n",
    "    return pmf_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`PmfOfWaitTime` makes a meta-Pmf that maps from each uniform\n",
    "distribution to its probability. Then it uses `MakeMixture`,\n",
    "which we saw in Section [mixture], to compute the mixture.\n",
    "\n",
    "`PmfOfWaitTime` also uses `MakeUniformPmf`,\n",
    "defined here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeUniformPmf(low, high):\n",
    "    pmf = thinkbayes.Pmf()\n",
    "    for x in MakeRange(low=low, high=high):\n",
    "        pmf.Set(x, 1)\n",
    "    pmf.Normalize()\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`low` and `high` are the range of the uniform\n",
    "distribution, (both ends included). Finally, `MakeUniformPmf`\n",
    "uses `MakeRange`, defined here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeRange(low, high, skip=10):\n",
    "    return range(low, high+skip, skip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`MakeRange` defines a set of possible values for wait time\n",
    "(expressed in seconds). By default it divides the range into 10 second\n",
    "intervals.\n",
    "\n",
    "To encapsulate the process of computing these distributions, I created a\n",
    "class called `WaitTimeCalculator`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class WaitTimeCalculator(object):\n",
    "\n",
    "    def __init__(self, pmf_z):\n",
    "        self.pmf_z = pmf_z\n",
    "        self.pmf_zb = BiasPmf(pmf)\n",
    "\n",
    "        self.pmf_y = self.PmfOfWaitTime(self.pmf_zb)\n",
    "        self.pmf_x = self.pmf_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameter, `pmf_z`, is the unbiased distribution of `z`.\n",
    "`pmf_zb` is the biased distribution of gap time, as seen by passengers.\n",
    "\n",
    "`pmf_y` is the distribution of wait time. `pmf_x` is the distribution of\n",
    "elapsed time, which is the same as the distribution of wait time. To see\n",
    "why, remember that for a particular value of `zp`, the\n",
    "distribution of `y` is uniform from 0 to `zp`.\n",
    "Also"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = zp - y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the distribution of `x` is also uniform from 0 to\n",
    "`zp`.\n",
    "\n",
    "Figure [fig.redline2] shows the distribution of `z`,\n",
    "`zb`, and `y` based on the data I collected from\n",
    "the Red Line web site.\n",
    "\n",
    "To present these distributions, I am switching from Pmfs to Cdfs. Most\n",
    "people are more familiar with Pmfs, but I think Cdfs are easier to\n",
    "interpret, once you get used to them. And if you want to plot several\n",
    "distributions on the same axes, Cdfs are the way to go.\n",
    "\n",
    "The mean of `z` is 7.8 minutes. The mean of `zb`\n",
    "is 8.8 minutes, about 13% higher. The mean of `y` is 4.4,\n",
    "half the mean of `zb`.\n",
    "\n",
    "As an aside, the Red Line schedule reports that trains run every 9\n",
    "minutes during peak times. This is close to the average of\n",
    "`zb`, but higher than the average of `z`. I\n",
    "exchanged email with a representative of the MBTA, who confirmed that\n",
    "the reported time between trains is deliberately conservative in order\n",
    "to account for variability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predicting wait times\n",
    "\n",
    "![Prior and posterior of `x` and predicted `y`.\n",
    "](figs/redline3.pdf)\n",
    "\n",
    "[fig.redline3]\n",
    "\n",
    "Let’s get back to the motivating question: suppose that when I arrive at\n",
    "the platform I see 10 people waiting. How long should I expect to wait\n",
    "until the next train arrives?\n",
    "\n",
    "As always, let’s start with the easiest version of the problem and work\n",
    "our way up. Suppose we are given the actual distribution of\n",
    "`z`, and we know that the passenger arrival rate, $\\lambda$,\n",
    "is 2 passengers per minute.\n",
    "\n",
    "In that case we can:\n",
    "\n",
    "1.  Use the distribution of `z` to compute the prior\n",
    "    distribution of `zp`, the time between trains as seen by\n",
    "    a passenger.\n",
    "\n",
    "2.  Then we can use the number of passengers to estimate the\n",
    "    distribution of `x`, the elapsed time since the last\n",
    "    train.\n",
    "\n",
    "3.  Finally, we use the relation `y = zp - x` to get the\n",
    "    distribution of `y`.\n",
    "\n",
    "The first step is to create a `WaitTimeCalculator` that\n",
    "encapsulates the distributions of `zp`, `x`, and\n",
    "`y`, prior to taking into account the number of passengers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "wtc = WaitTimeCalculator(pmf_z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pmf_z` is the given distribution of gap times.\n",
    "\n",
    "The next step is to make an `ElapsedTimeEstimator` (defined\n",
    "below), which encapsulates the posterior distribution of `x`\n",
    "and the predictive distribution of `y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ete = ElapsedTimeEstimator(wtc,\n",
    "                           lam=2.0/60,\n",
    "                           num_passengers=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters are the `WaitTimeCalculator`, the passenger\n",
    "arrival rate, `lam` (expressed in passengers per second), and\n",
    "the observed number of passengers, let’s say 15.\n",
    "\n",
    "Here is the definition of `ElapsedTimeEstimator`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class ElapsedTimeEstimator(object):\n",
    "\n",
    "    def __init__(self, wtc, lam, num_passengers):\n",
    "        self.prior_x = Elapsed(wtc.pmf_x)\n",
    "\n",
    "        self.post_x = self.prior_x.Copy()\n",
    "        self.post_x.Update((lam, num_passengers))\n",
    "\n",
    "        self.pmf_y = PredictWaitTime(wtc.pmf_zb, self.post_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`prior_x` and `posterior_x` are the prior and posterior distributions of\n",
    "elapsed time. `pmf_y` is the predictive distribution of wait time.\n",
    "\n",
    "`ElapsedTimeEstimator` uses `Elapsed` and\n",
    "`PredictWaitTime`, defined below.\n",
    "\n",
    "`Elapsed` is a Suite that represents the hypothetical\n",
    "distribution of `x`. The prior distribution of `x`\n",
    "comes straight from the `WaitTimeCalculator`. Then we use the\n",
    "data, which consists of the arrival rate, `lam`, and the\n",
    "number of passengers on the platform, to compute the posterior\n",
    "distribution.\n",
    "\n",
    "Here’s the definition of `Elapsed`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Elapsed(thinkbayes.Suite):\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        x = hypo\n",
    "        lam, k = data\n",
    "        like = thinkbayes.EvalPoissonPmf(k, lam * x)\n",
    "        return like"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, `Likelihood` takes a hypothesis and data, and\n",
    "computes the likelihood of the data under the hypothesis. In this case\n",
    "`hypo` is the elapsed time since the last train and\n",
    "`data` is a tuple of `lam` and the number of\n",
    "passengers.\n",
    "\n",
    "The likelihood of the data is the probability of getting `k`\n",
    "arrivals in `x` time, given arrival rate `lam`. We\n",
    "compute that using the PMF of the Poisson distribution.\n",
    "\n",
    "Finally, here’s the definition of `PredictWaitTime`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def PredictWaitTime(pmf_zb, pmf_x):\n",
    "    pmf_y = pmf_zb - pmf_x\n",
    "    RemoveNegatives(pmf_y)\n",
    "    return pmf_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pmf_zb` is the distribution of gaps between trains; `pmf_x` is the\n",
    "distribution of elapsed time, based on the observed number of\n",
    "passengers. Since `y = zb - x`, we can compute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pmf_y = pmf_zb - pmf_x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The subtraction operator invokes `Pmf.__sub__`, which enumerates all\n",
    "pairs of `zb` and `x`, computes the differences,\n",
    "and adds the results to `pmf_y`.\n",
    "\n",
    "The resulting Pmf includes some negative values, which we know are\n",
    "impossible. For example, if you arrive during a gap of 5 minutes, you\n",
    "can’t wait more than 5 minutes. `RemoveNegatives` removes the\n",
    "impossible values from the distribution and renormalizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RemoveNegatives(pmf):\n",
    "    for val in pmf.Values():\n",
    "        if val < 0:\n",
    "            pmf.Remove(val)\n",
    "    pmf.Normalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure [fig.redline3] shows the results. The prior distribution of\n",
    "`x` is the same as the distribution of `y` in\n",
    "Figure [fig.redline2]. The posterior distribution of `x`\n",
    "shows that, after seeing 15 passengers on the platform, we believe that\n",
    "the time since the last train is probably 5-10 minutes. The predictive\n",
    "distribution of `y` indicates that we expect the next train\n",
    "in less than 5 minutes, with about 80% confidence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating the arrival rate\n",
    "\n",
    "![Prior and posterior distributions of `lam` based on five\n",
    "days of passenger data. ](figs/redline1.pdf)\n",
    "\n",
    "[fig.redline1]\n",
    "\n",
    "The analysis so far has been based on the assumption that we know (1)\n",
    "the distribution of gaps and (2) the passenger arrival rate. Now we are\n",
    "ready to relax the second assumption.\n",
    "\n",
    "Suppose that you just moved to Boston, so you don’t know much about the\n",
    "passenger arrival rate on the Red Line. After a few days of commuting,\n",
    "you could make a guess, at least qualitatively. With a little more\n",
    "effort, you could estimate $\\lambda$ quantitatively.\n",
    "\n",
    "Each day when you arrive at the platform, you should note the time and\n",
    "the number of passengers waiting (if the platform is too big, you could\n",
    "choose a sample area). Then you should record your wait time and the\n",
    "number of new arrivals while you are waiting.\n",
    "\n",
    "After five days, you might have data like this:\n",
    "\n",
    "    k1      y     k2\n",
    "    --     ---    --\n",
    "    17     4.6     9\n",
    "    22     1.0     0\n",
    "    23     1.4     4\n",
    "    18     5.4    12\n",
    "    4      5.8    11\n",
    "\n",
    "where `k1` is the number of passengers waiting when you\n",
    "arrive, `y` is your wait time in minutes, and `k2`\n",
    "is the number of passengers who arrive while you are waiting.\n",
    "\n",
    "Over the course of one week, you waited 18 minutes and saw 36 passengers\n",
    "arrive, so you would estimate that the arrival rate is 2 passengers per\n",
    "minute. For practical purposes that estimate is good enough, but for the\n",
    "sake of completeness I will compute a posterior distribution for\n",
    "$\\lambda$ and show how to use that distribution in the rest of the\n",
    "analysis.\n",
    "\n",
    "`ArrivalRate` is a `Suite` that represents\n",
    "hypotheses about $\\lambda$. As always, `Likelihood` takes a\n",
    "hypothesis and data, and computes the likelihood of the data under the\n",
    "hypothesis.\n",
    "\n",
    "In this case the hypothesis is a value of $\\lambda$. The data is a pair,\n",
    "`y, k`, where `y` is a wait time and\n",
    "`k` is the number of passengers that arrived."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class ArrivalRate(thinkbayes.Suite):\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        lam = hypo\n",
    "        y, k = data\n",
    "        like = thinkbayes.EvalPoissonPmf(k, lam * y)\n",
    "        return like"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This `Likelihood` might look familiar; it is almost identical\n",
    "to `Elapsed.Likelihood` in Section [elapsed]. The difference\n",
    "is that in `Elapsed.Likelihood` the hypothesis is\n",
    "`x`, the elapsed time; in `ArrivalRate.Likelihood`\n",
    "the hypothesis is `lam`, the arrival rate. But in both cases\n",
    "the likelihood is the probability of seeing `k` arrivals in\n",
    "some period of time, given `lam`.\n",
    "\n",
    "`ArrivalRateEstimator` encapsulates the process of estimating\n",
    "$\\lambda$. The parameter, `passenger_data`, is a list of `k1, y,\n",
    "k2` tuples, as in the table above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class ArrivalRateEstimator(object):\n",
    "\n",
    "    def __init__(self, passenger_data):\n",
    "        low, high = 0, 5\n",
    "        n = 51\n",
    "        hypos = numpy.linspace(low, high, n) / 60\n",
    "\n",
    "        self.prior_lam = ArrivalRate(hypos)\n",
    "\n",
    "        self.post_lam = self.prior_lam.Copy()\n",
    "        for k1, y, k2 in passenger_data:\n",
    "            self.post_lam.Update((y, k2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`__init__` builds `hypos`, which is a sequence of\n",
    "hypothetical values for `lam`, then builds the prior\n",
    "distribution, `prior_lam`. The `for` loop updates the prior\n",
    "with data, yielding the posterior distribution, `post_lam`.\n",
    "\n",
    "Figure [fig.redline1] shows the prior and posterior distributions. As\n",
    "expected, the mean and median of the posterior are near the observed\n",
    "rate, 2 passengers per minute. But the spread of the posterior\n",
    "distribution captures our uncertainty about $\\lambda$ based on a small\n",
    "sample."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Incorporating uncertainty\n",
    "\n",
    "![Predictive distributions of `y` for possible values of\n",
    "`lam`. ](figs/redline4.pdf)\n",
    "\n",
    "[fig.redline4]\n",
    "\n",
    "Whenever there is uncertainty about one of the inputs to an analysis, we\n",
    "can take it into account by a process like this:\n",
    "\n",
    "1.  Implement the analysis based on a deterministic value of the\n",
    "    uncertain parameter (in this case $\\lambda$).\n",
    "\n",
    "2.  Compute the distribution of the uncertain parameter.\n",
    "\n",
    "3.  Run the analysis for each value of the parameter, and generate a set\n",
    "    of predictive distributions.\n",
    "\n",
    "4.  Compute a mixture of the predictive distributions, using the weights\n",
    "    from the distribution of the parameter.\n",
    "\n",
    "We have already done steps (1) and (2). I wrote a class called\n",
    "`WaitMixtureEstimator` to handle steps (3) and (4)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "class WaitMixtureEstimator(object):\n",
    "\n",
    "    def __init__(self, wtc, are, num_passengers=15):\n",
    "        self.metapmf = thinkbayes.Pmf()\n",
    "\n",
    "        for lam, prob in sorted(are.post_lam.Items()):\n",
    "            ete = ElapsedTimeEstimator(wtc, lam, num_passengers)\n",
    "            self.metapmf.Set(ete.pmf_y, prob)\n",
    "\n",
    "        self.mixture = thinkbayes.MakeMixture(self.metapmf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`wtc` is the `WaitTimeCalculator` that contains\n",
    "the distribution of `zb`. `are` is the\n",
    "`ArrivalTimeEstimator` that contains the distribution of\n",
    "`lam`.\n",
    "\n",
    "The first line makes a meta-Pmf that maps from each possible\n",
    "distribution of `y` to its probability. For each value of\n",
    "`lam`, we use `ElapsedTimeEstimator` to compute\n",
    "the corresponding distribution of `y` and store it in the\n",
    "Meta-Pmf. Then we use `MakeMixture` to compute the mixture.\n",
    "\n",
    "Figure [fig.redline4] shows the results. The shaded lines in the\n",
    "background are the distributions of `y` for each value of\n",
    "`lam`, with line thickness that represents likelihood. The\n",
    "dark line is the mixture of these distributions.\n",
    "\n",
    "In this case we could get a very similar result using a single point\n",
    "estimate of `lam`. So it was not necessary, for practical\n",
    "purposes, to include the uncertainty of the estimate.\n",
    "\n",
    "In general, it is important to include variability if the system\n",
    "response is non-linear; that is, if small changes in the input can cause\n",
    "big changes in the output. In this case, posterior variability in\n",
    "`lam` is small and the system response is approximately\n",
    "linear for small perturbations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decision analysis \n",
    "\n",
    "![Probability that wait time exceeds 15 minutes as a function of the\n",
    "number of passengers on the platform. ](figs/redline5.pdf)\n",
    "\n",
    "[fig.redline5]\n",
    "\n",
    "At this point we can use the number of passengers on the platform to\n",
    "predict the distribution of wait times. Now let’s get to the second part\n",
    "of the question: when should I stop waiting for the train and go catch a\n",
    "taxi?\n",
    "\n",
    "Remember that in the original scenario, I am trying to get to South\n",
    "Station to catch the commuter rail. Suppose I leave the office with\n",
    "enough time that I can wait 15 minutes and still make my connection at\n",
    "South Station.\n",
    "\n",
    "In that case I would like to know the probability that `y`\n",
    "exceeds 15 minutes as a function of `num_passengers`. It is easy enough\n",
    "to use the analysis from Section [elapsed] and run it for a range of\n",
    "`num_passengers`.\n",
    "\n",
    "But there’s a problem. The analysis is sensitive to the frequency of\n",
    "long delays, and because long delays are rare, it is hard estimate their\n",
    "frequency.\n",
    "\n",
    "I only have data from one week, and the longest delay I observed was 15\n",
    "minutes. So I can’t estimate the frequency of longer delays accurately.\n",
    "\n",
    "However, I can use previous observations to make at least a coarse\n",
    "estimate. When I commuted by Red Line for a year, I saw three long\n",
    "delays caused by a signaling problem, a power outage, and “police\n",
    "activity” at another stop. So I estimate that there are about 3 major\n",
    "delays per year.\n",
    "\n",
    "But remember that my observations are biased. I am more likely to\n",
    "observe long delays because they affect a large number of passengers. So\n",
    "we should treat my observations as a sample of `zb` rather\n",
    "than `z`. Here’s how we can do that.\n",
    "\n",
    "During my year of commuting, I took the Red Line home about 220 times.\n",
    "So I take the observed gap times, `gap_times`, generate a sample of 220\n",
    "gaps, and compute their Pmf:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import thinkbayes\n",
    "\n",
    "n = 220\n",
    "cdf_z = thinkbayes.MakeCdfFromList(gap_times)\n",
    "sample_z = cdf_z.Sample(n)\n",
    "pmf_z = thinkbayes.MakePmfFromList(sample_z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next I bias `pmf_z` to get the distribution of `zb`, draw a\n",
    "sample, and then add in delays of 30, 40, and 50 minutes (expressed in\n",
    "seconds):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cdf_zp = BiasPmf(pmf_z).MakeCdf()\n",
    "sample_zb = cdf_zp.Sample(n) + [1800, 2400, 3000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Cdf.Sample` is more efficient than `Pmf.Sample`,\n",
    "so it is usually faster to convert a Pmf to a Cdf before sampling.\n",
    "\n",
    "Next I use the sample of `zb` to estimate a Pdf using KDE,\n",
    "and then convert the Pdf to a Pmf:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pdf_zb = thinkbayes.EstimatedPdf(sample_zb)\n",
    "xs = MakeRange(low=60)\n",
    "pmf_zb = pdf_zb.MakePmf(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally I unbias the distribution of `zb` to get the\n",
    "distribution of `z`, which I use to create the\n",
    "`WaitTimeCalculator`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pmf_z = UnbiasPmf(pmf_zb)\n",
    "wtc = WaitTimeCalculator(pmf_z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This process is complicated, but all of the steps are operations we have\n",
    "seen before. Now we are ready to compute the probability of a long wait."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ProbLongWait(num_passengers, minutes):\n",
    "    ete = ElapsedTimeEstimator(wtc, lam, num_passengers)\n",
    "    cdf_y = ete.pmf_y.MakeCdf()\n",
    "    prob = 1 - cdf_y.Prob(minutes * 60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the number of passengers on the platform,\n",
    "`ProbLongWait` makes an `ElapsedTimeEstimator`,\n",
    "extracts the distribution of wait time, and computes the probability\n",
    "that wait time exceeds `minutes`.\n",
    "\n",
    "Figure [fig.redline5] shows the result. When the number of passengers is\n",
    "less than 20, we infer that the system is operating normally, so the\n",
    "probability of a long delay is small. If there are 30 passengers, we\n",
    "estimate that it has been 15 minutes since the last train; that’s longer\n",
    "than a normal delay, so we infer that something is wrong and expect\n",
    "longer delays.\n",
    "\n",
    "If we are willing to accept a 10% chance of missing the connection at\n",
    "South Station, we should stay and wait as long as there are fewer than\n",
    "30 passengers, and take a taxi if there are more.\n",
    "\n",
    "Or, to take this analysis one step further, we could quantify the cost\n",
    "of missing the connection and the cost of taking a taxi, then choose the\n",
    "threshold that minimizes expected cost."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "The analysis so far has been based on the assumption that the arrival\n",
    "rate of passengers is the same every day. For a commuter train during\n",
    "rush hour, that might not be a bad assumption, but there are some\n",
    "obvious exceptions. For example, if there is a special event nearby, a\n",
    "large number of people might arrive at the same time. In that case, the\n",
    "estimate of `lam` would be too low, so the estimates of\n",
    "`x` and `y` would be too high.\n",
    "\n",
    "If special events are as common as major delays, it would be important\n",
    "to include them in the model. We could do that by extending the\n",
    "distribution of `lam` to include occasional large values.\n",
    "\n",
    "We started with the assumption that we know distribution of\n",
    "`z`. As an alternative, a passenger could estimate\n",
    "`z`, but it would not be easy. As a passenger, you only\n",
    "observe only your own wait time, `y`. Unless you skip the\n",
    "first train and wait for the second, you don’t observe the gap between\n",
    "trains, `z`.\n",
    "\n",
    "However, we could make some inferences about `zb`. If we note\n",
    "the number of passengers waiting when we arrive, we can estimate the\n",
    "elapsed time since the last train, `x`. Then we observe\n",
    "`y`. If we add the posterior distribution of `x`\n",
    "to the observed `y`, we get a distribution that represents\n",
    "our posterior belief about the observed value of `zb`.\n",
    "\n",
    "We can use this distribution to update our beliefs about the\n",
    "distribution of `zb`. Finally, we can compute the inverse of\n",
    "`BiasPmf` to get from the distribution of `zb` to\n",
    "the distribution of `z`.\n",
    "\n",
    "I leave this analysis as an exercise for the reader. One suggestion: you\n",
    "should read Chapter [species] first. You can find the outline of a\n",
    "solution in <http://thinkbayes.com/redline.py>. For more information see\n",
    "Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "\n",
    "This exercise is from MacKay, *Information Theory, Inference, and\n",
    "Learning Algorithms*:\n",
    "\n",
    "> Unstable particles are emitted from a source and decay at a distance\n",
    "> $x$, a real number that has an exponential probability distribution\n",
    "> with [parameter] $\\lambda$. Decay events can only be observed if they\n",
    "> occur in a window extending from $x=1$ cm to $x=20$ cm. $N$ decays are\n",
    "> observed at locations $\\{ 1.5, 2, 3, 4, 5, 12 \\}$ cm. What is the\n",
    "> posterior distribution of $\\lambda$?\n",
    "\n",
    "You can download a solution to this exercise from\n",
    "<http://thinkbayes.com/decay.py>."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
